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Abstract 

The study of planktonic ecosystems is important as they make up the 
bottom trophic levels of aquatic food webs. We study a closed Nutrient- 
Phytoplankton-Zooplankton (NPZ) model that includes size structure in 
the juvenile zooplankton. The closed nature of the system allows the 
formulation of a conservation law of biomass that governs the system. 
The model consists of a system of nonlinear ordinary differential equation 
coupled to a partial differential equation. We are able to transform this 
system into a system of delay differential equations where the delay is of 
threshold type and is state-dependent. The system of delay differential 
equations can be further transformed into one with fixed delay. Using 
the different forms of the model we perform a qualitative analysis of the 
solutions, which includes studying existence and uniqueness, positivity 
and boundedness, local and global stability, and conditions for extinction. 
Key parameters that are explored are the total biomass in the system 
and the maturity level at which the juvenile zooplankton reach maturity. 
Numerical simulations are also performed to verify our analytical results. 


1 Introduction 

Nutrient-Phytoplankton-Zooplankton (NPZ) models are used to describe the 
bottom two trophic levels of an aquatic ecosystem. As is the case with many 
ecological models, they range from very simple to very complex. Simple models, 
such as the Lotka-Volterra system ns, are beneficial in that one can obtain 
analytical results more easily, but often suffer from lack of realism. On the other 
hand, complex models may theoretically represent a more accurate description 
of reality, but may be difficult or impossible to understand in any general way, 
and may be useless without precise and accurate parameter values. A lot of 
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structure can be given to an NPZ model, such as size-dependence and spatial 
dependence, and this structure leads to rich modelling possibilities such as size- 
dependent parameters and coupling with external factors like fluid dynamics 
or higher predation. To strike a balance, it may be useful to start with a 
simple model that focuses on one factor affecting the ecosystem, and use it to 
obtain analytical results and study general trends. For an introduction in the 
construction of NPZ models, see [6j. 

We will focus on the role of maturity in the juvenile zooplankton population 
within an ecosystem while keeping all other factors as simple as possible. In 
this vain, a simple NPZ model is coupled with a standard linear first-order PDE 
that describes the spectrum of the juvenile zooplankton population as a func¬ 
tion of time and maturity. It is a known result that this type of PDE is closely 
related to delay equations. Examples of population models where time delay is 
a consequence of age structure can be found in We will be 

considering size structure where the rate of growth of the juvenile zooplankton 
is permitted to depend on the concentration of the phytoplankton. With this 
assumption, the delay in the related delay equation is not explicitly defined, but 
rather defined implicitly through a threshold-type condition. Consequently, sys¬ 
tems with this type of delay are known as threshold delay differential equations. 
They have been studied in the case of a single population in [2D] and for an 
insect species in 13, for example. These models operate under the assumption 
that maturation occurs when an immature individual accumulates enough of 
some quantity, such as size or weight. Other applications include red blood cell 
production [15] and immunology [251 . 

Our model is formulated in such a way that biomass is conserved. In other 
words, we are considering a closed ecosystem with no mass being added to or 
subtracted from the system. NPZ models with this property have been studied, 
for example, by [7] and [26) . A typical property of these systems is that the 
amount of biomass, which is determined by initial conditions, plays a crucial 
role in determining the types of dynamics that can and do occur. Most notably, 
an insufficient amount of biomass leads to the extinction of plankton; either the 
zooplankton only, or both the phytoplankton and zooplankton. 


2 Structured Model 

We consider a Nutrient-Phytoplankton-Zooplankton (NPZ) model in which the 
zooplankton population is split into a mature class and a juvenile class. We 
will consider the juvenile class as a function of time and maturity, denoted 
Maturity is considered to be an abstract quantity, of which the juvenile 
zooplankton much accumulate enough in order to enter adulthood. The total 
concentration of immature zooplankton with maturity levels between si and S 2 
at time t is then J s * 2 p(t, s ) ds. We restrict s to the interval [0, m], where m is 
the required level of maturity for adulthood. The nutrient, phytoplankton, and 
mature zooplankton variables will therefore only depend on time. 

Since the zooplankton feed on phytoplankton, we will assume that the growth 
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rate of the juvenile zooplankton depends on the quantity of phytoplankton 
present in the system. We will denote this dependence as R(P(t)). 

The following equations model the ecosystem and we will refer to it as the 
PDE1 model. 


(la) 

(lb) 

(lc) 

(l d) 

(le) 


= - pP(t)f(N(t )) + A P{t) + 5Z(t) + (1 - 'y)gZ(t)h(P(t)) 

r-m 

+ / 5 0 p(t,s)ds, 

Jo 

^ = pP(t)f(N(t)) - XP{t) - gZ(t)h(P(t)) : 

= ^(- p (*))p(*> m ) - sz(t), 

s) + R{P(t))^(t, s) = - 8 0 p{t , s), 

R(P(t))p(t, o) =7 gZ(t)h(P(t)). 


Appropriate initial conditions for these equations are 

(2) N(0) = N 0 , P(0)=P o , Z(0) = Z 0 , p(0 ,s) = p 0 (s), 

where Nq , Pq , and Zq are nonnegative real numbers and po is a nonnegative 
function on the interval [0,m]. 

A few basic ecological processes govern the system. Phytoplankton, P, up¬ 
take nutrient, N , at a rate that is proportional to the quantity of phytoplankton 
and a function of the dissolved nutrient, /(IV). Zooplankton, Z, graze on phy¬ 
toplankton at a rate that is proportional to the quantity of zooplankton and a 
function of phytoplankton, h(P). This is marked by a grazing efficiency factor, 
7 £ (0,1]. Phytoplankton mortality is proportional to the amount of phyto¬ 
plankton present at the time. Zooplankton mortality is similarly proportional 
to the amount of zooplankton in the system at that time. Nutrient recycling 
transforms dead biomass and faecal matter back to the dissolved nutrient vari¬ 
able. Equation (lldl) is a standard transport equation [T5] with a decay term 
due to natural mortality. The boundary condition, given in equation m, says 
that the birth rate of the immature zooplankton is equal to some fraction of the 
biomass obtained by the zooplankton through grazing. 

Some assumptions are made on the functional responses. For the phyto¬ 
plankton nutrient uptake response, we assume / £ C 1 and 

(3) /(0) = 0, /'(TV) > 0, lim /(TV) = 1. 

iV-j-oo 

Similarly, for the zooplankton functional response for grazing on phytoplankton, 
we assume h £ C 1 and 


(4) h{ 0) = 0, h\P) > 0, lim h{P) = 1. 

P—YOO 

These assumptions on h encompass both Type II (concave down) and Type III 
(sigmoidal shaped) responses, as described in . 
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We will also assume throughout that all parameter values are positive (al¬ 
though we will allow <5 q to be zero) and that 


M > A, 7 g> 8, 

so that / _1 ( X/p) and h~ l {5/{^g)) exist. 

It will be assumed that R £ C 1 and that it satisfies the following properties. 


(5) R(P) > 0, R'(P) > 0, R'{ 0) > 0 if i?(0) = 0, lim R(P) = R x < oo. 

P—>oo 


These assumptions can be justified by laboratory experiments in [14] . where 
they measured the development rate of the zooplankter Daphnia as a function 
of F, the amount of available food. They found that the development rate was 
proportional to F/(F + Fhaif )■ 

The assumptions on h in equation (J3J and the third assumption on R in 
equation (JSJ) imply that 


( 6 ) 


lim W 
p-> o+ R(P) 


< oo. 


This puts a bound on p(f, 0) as the phytoplankton population approaches zero, 
as seen in equation ©• 

As previously stated, there is no biomass lost from system 0. We can 
confirm this by noting that PDE1 model 0 satisfies the following conservation 
law: 

pm 

(7) N t = N(t) + P(t) T Z(t) + / p{t, s) ds, 

Jo 

where Nt is the total biomass in the system. Given this conservation law, the 
total amount of juvenile zooplankton can be determined from the total biomass 
and the quantity of the dissolved nutrient, phytoplankton, and zooplankton. 
Making the appropriate substitution of equation 0 into 0 yields the following 
system: 


(8a) 

(8b) 

(8c) 

(8d) 

(8e) 


^ = - pP(t)f(N(t )) + A P{t) + SZ(t) + (1 - 7 )gZ{t)h(P(t)) 
+ 6 0 (N T -N(t)-P(t)-Z(t)), 

^ = pP(t)f(N(t)) - A P{t) - gZ(t)h(P(t)), 

= RiPifypiti m ) - sz (t), 

s) + R{P(t))^-(t, s) = -6 0 p(t, s), 

R(P(t))p(t, o) =7 gZ(t)h(P(t)). 


We will refer to this system as the PDE2 model. While the substitution is 
straightforward, it is important to note that we have now fixed Nt- Thus, in 
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order for a solution of system © to be the same as a solution to system © , we 
must choose initial conditions in @ that satisfy 

pm 

(9) Nt = N 0 + Pq + Z 0 + / p 0 (s)ds. 

J o 

In other words, the total biomass, Nt, is set by the initial conditions for the 
PDE1 model ©, so if we treat it as a parameter in the PDE2 model ©, then 
we have to restrict the possible initial conditions by the relation (|TJ|) . If we fix 
Nt in system © and choose initial conditions that do not satisfy ©, then the 
resulting solution is not a solution to PDE1. Indeed, if we define 

pm 

(10) A(t) = — Nt H - N(t ) + P(t) + Z(t) / p(t , s) ds, 

Jo 

it can be verified that a solution to © satisfies 

(11) j t A (t) = -5 0 A(t). 

This implies that if equation © is not satisfied (A(0) ^ 0), then equation © 
is never satisfied (A (t) ^ 0). However, a solution for PDE1 always satisfies 
equation © (A (t) = 0). We can conclude, however, that a solution to system 
PDE2 will asymptotically approach a solution to PDE1, assuming that both 
solutions exist for all time. Due to the simpler nature of the PDE2 model, we 
will study it instead of the original system, and assume that initial conditions 
satisfy equation ©• 

[20] studies a model in a similar form to the PDE2 model ©, although the 
system is for a single species where the mature population is coupled with its 
own immature class. The growth rate of the immature class depended on the 
size of the mature population. m reduced the equations so that they became 
a scalar threshold delay equation. This type of equation has a state-dependent 
delay due to the fact that a new-born individual must first reach maturity 
before entering into the adult population. This delay is defined implicitly via 
an integral which represents the condition for reaching maturity. [20] then uses 
a clever change of variables to transform the system to a functional differential 
equation where the delay no longer depends on the state. Due to the similar 
form of © to Smith’s model, we will use the same methods to reduce our model 
to a threshold delay equation and then to a delay differential equation with a 
bounded and state-independent delay. 


3 Reduction to a Threshold Delay Equation 


As in [20], we can solve (I8dl) and (l8cl) for p as a function of P and Z using the 
method of characteristics, and then substitute this solution into (IScl) . It can be 
verified that the solution for p is 


(12) p(t,s) 


e Sot p 0 (s - f* R(P(u)) duj 

r -SnT(s.P t ) jgZ(t-T(a,Pt))h(P(t-T(s,Pt))) 
e R(P(t-r(s,P t ))) 


for (t, s ) G Si 
for (f, s) € S 2 
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where Pt(u) = P{t + u), and r(s, Pt) is defined implicitly by 

(13) f R(P(t + u)) du = s. 

J-T(s,P t ) 

Biologically, it is the time it takes for the immature zooplankton to grow to size s 
as a function of the history of phytoplankton. For example, if the phytoplankton 
population is at an equilibrium value Pit) = P* for all t, then r(s, P t ) = s p . 

Substituting p(t,m), defined by equation (fl2]l . into <|8cl> gives a nonau- 
tonomous threshold delay equation. Alternatively, we can change our initial 
conditions from fl2J) to 

(14) N(to + t) = 4>i{t), P(to + t) = fait), Z(to + t) = <f>z(t ), 

for t £ [—to, 0], with t(to, 0 2 ) = to- The initial conditions are shifted in a 
non-standard way so they correspond with the original PDE1 model £[]). The 
time to is when all of the immature zooplankton initially present at t = 0 have 
reached maturity. With these initial conditions, we can restrict p(t,s) to S 2 - 
Then we have the following autonomous system 

^ = - pP(t)f(N(t )) + A Pit) + SZ(t) + (1 - 7 )gZ(t)h(P(t)) 

(15a) + S 0 (N t - N(t) - Pit) — Z(t)), 

(15b) 

^ =»P(t)f(N(t)) - AP(t) - gZ(t)h(P(t)), 

(15c) 

<m. =R{m)e -s 0 r M) IdZit - P t ))) _ 

dt y R(P(t-T(m,P t ))) 

for t > t 0 , where the function r is given by equation m- We will refer to this 
system as the TDE model, as systems in this form are known as threshold delay 
equations. 

Substituting p(t,s) for (t, s) £ S 2 into the definition of A in equation (flOll . 
we find that 


A(t) = — Nt + N(t) + Pit) + Z(t) 


(16) 


f 


„-s 0 T( a ,p t ) 7 gZ (t - r(s, P t ))h{P(t - r(s, P t ))) 
R(P(t-Ti Sl P t ))) 


ds. 


It can be verified that solutions to system (fl5l) still satisfy equation (fill) with 
A defined by equation (fTBl) . In particular, if the initial conditions (fill) satisfy 


Np =^i (0) + 02(0) + 03(0) 


(17) 


—<5 o t(s, 0 2 ) 7g03(-T(a, 0 2 ))h(0 2 (—t(s, 0 2 ))) 

R(M-T(s,<h))) 
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then we have 


(18) 


N t =N(t) + P(t) + Z(t) 



-SoT(s,Pt) 7 9 z (t ~ T 0, P t ))h(P(t - t(s, P t ))) 
R(P(t — t(s, P t ))) 


ds : 


for t > tg. Since initial conditions satisfying equation m are the only ones that 
relate to the original system ©, we will restrict our attention to this particular 
choice. For later use we define the set 


D Nt 3 ) t e C([-t o ,0],R 3 ) : cj)i(t ) > 0,i = 1,2,3, t e [-i 0 ,0], 

r(m, 4 > 2 ) = to, and (flTl) is satisfied}. 

This is the set of initial conditions that we will consider. That is, we will not 
study cases where N, P, or Z have initial values equal to zero as we are interested 
in studying the dynamics when all three trophic levels are initially present. 

We note that system m is a differential equation with locally bounded 
delay. Since we allow for the case where R( 0) = 0, it is possible that the 
delay becomes arbitrarily large if the phytoplankton population approaches zero. 
However, for any given state, the delay is bounded for all states within its 
neighbourhood. Results for differential equations with locally bounded delay can 
be found in [24]. However, our system is in the form that it my be transformed 
into a system with state-independent delay as in [19j [20]. We will therefore use 
this approach. 


4 Transformation to a State-Independent Delay 
Differential Equation 

It is possible to remove the state-dependence on the delay through a clever 
change in variable, as given in uansi. The transformation is a state-dependent 
change of the time variable: 


(19a) 

r r R{P{u)) , 

t= —- du, 

Jo R* 

(19b) 

SO 

II 

<So' 

(19c) 

Hi) = P(t), 

(19d) 

II 

cr-s 


for t > 0, where R* is a typical value of R{P). With this transformation 
t — t(s , P t ) becomes a constant, since for t > s/R* we have 



r r(p(u)) du _ r r(p(u)) du = /■ t - T(s ’ Pt) r(p(u)) 
Jo R* Jt-T(s,p t ) R* Jo R* 
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from the definition of r in equation (fl3l) . Then P(t — s/R*) = P(t — r(s, Pt)) 
and Z(t — s/R*) = Z(t — r(s,Pt)). That is, the delay no longer depends on 
the state. Furthermore, we can obtain an explicit equation for r. Define f : 
[0, m] x C[—m/R*,0] —> [0, oo) through 


( 20 ) 


t(s, Pf) = t(s, P t ) = du = 


R* 


t—T(s,Pt) Ji—s/R* R{P{r)) 


■ dr , 


where we have used the substitution 


R(P(y)) 

R* 


dv. 


In a similar way, we can invert the time transformation: 


( 21 ) 



/'* R* 

/ - 7 - dr. 

Jo R(P(r)) 


With this change of variables, the TDE model (fTSl) is transformed to the 
following delay differential equation: 


=1+ xPp ) + + C 1 - 7 )gZ(i)h(P(t)) 

dt R{P{t)) 

(22a) + S 0 (N t - N(t) - Pit) - Z(t))], 

(22b) 

= r<M) 1 >^(^ (*)) xP (*) - 9Z p )h( p dm 


(22c) 
dZjt ) 
dt 




R* 


R(P{t-T)) 


Z(t — T)h(P{t — T)) - 


R* 


Rim) 


5Z(t), 


for t >T, where T = m/R* and f (m, P/) is defined by equation (EOl) . We will 
refer to this system as the DDE model. The corresponding initial conditions 
are 


(23) N(T + t) = 4> 1 (t), P(T + t) = fo(t), Z(T + t)=fo(t), 


for t £ [—T, 0], where 0i,<fc,<^3 G C[—T, 0]. Again, solutions corresponding to 
the PDE1 model ([!]) should satisfy 

(24) Nt = </>i(0) + 02(0) + </> 3 (0) + [ dg _ 

Jo R{M~w)) 

For an initial condition that satisfies equation EH) , the solution to system (1221) 
satisfies 

(25) N t = N(t) + P(t) + Z(t) + f e - 5 °Rm ~ y ) h i Pp ~ w)) ds , 

Jo R{P{t - -^)) 
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for all t > T for which the solution exists. 

If we were to apply the time transformation in equation (119al) directly to 
the PDE2 model (|8j). and then proceed with the method of characteristics, 
we would arrive at the DDE model (1221) . In essence, the time transformation 
removes the dependency of the characteristic curves on P(t), and instead places 
this dependence on the solutions along the curves. In m and [22], he applies 
an analogous time transform to a system coupling a single adult population 
coupled with its corresponding maturity-structured juvenile population. 

Since the delay in the DDE model (1221) is bounded and state-independent, 
the standard theory for functional differential equations m may be applied. 


5 Existence and Uniqueness of Solutions 

In this section, we will develop conditions for there to be a unique solution to 
the initial value problem for the TDE model (fbH) - (fT5l) that exists for all time. 
This will be done by matching it to a corresponding initial value problem for 
the DDE model (1221) - (l23l) . First, define the open set 

(26) n = {(&,&, $,) T G C[-T, 0] : R(M0)) P 0 for 6 £ [-T,0]}. 

Standard results on functional differential equations with bounded delay (for 
example, see m ) give us that for (0i, <j> 2 ,03 ) T G Cl, the initial value problem 
(1221) - (1251l has a unique maximal solution. We can then deduce the following 
result. 

Proposition 5.1 If <j>i, fa, (f>3 are positive functions that satisfy equation (1241) . 
then the unique maximal solution to the initial value problem m-m, denoted 
as ( N,P,Z) T : [0,w) —> R 3 , satisfies N(t), P(t), Z(t) £ (0 ,Nt) for all t £ 
[0 ,w). 

Proof: Since 4>\,4>2,4>3 are positive functions, the result is clearly true for 

t £ [0, T], By the properties of R in equation ([5]) and the definition of the phase 
space, Cl in equation (1261) . we find that P[t) > 0 for all t £ [0, w). Suppose there 
exists fii> T such that Z(f3i) = 0 and Z{t) > 0 for t < /3i. Then, for t £ [T, /3i] 

we have that > —SZ(t), which implies that Z{i) >0 since Z(T) > 0. In 
particular, Z(f3i) > 0, which is a contradiction, so it must be true that Z(t) > 0 
for all t £ [0, w). Equation (1251) then implies that N(t)+P(t)+Z(t) < Nt, which 
implies that N is increasing at sufficiently small and positive values of N(t). 
Therefore, N(t) > 0 for all t £ [0,w). From equation (1251) . N(i),P(t), Z{i) > 0 
implies also that N(t), P(t), Z(t) < Nt- □ 

It follows from this Proposition that all solutions are bounded. From the 
definition of the phase space, Cl in equation (l26l) , standard results for functional 
differential equations m imply that each solution to the initial value problem 
(l22l) - (l2dl) either exists for all time or satisfies 

lim R(P(i)) = 0. 

t—yw 
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That is, solutions either exist for all time or approach dll. 

Proposition 5.2 If (<t>i,4 > 2, < t > 3) T E Pat t , then the solution to the initial value 
problem ©-(usd exists for all time and satisfies N(t),P(t),Z(t) £ (0,Nt) for 
all t £ [0, oo). 

Proof: Consider the corresponding initial value problem (l22l) - (l23l) with (f>i(t) = 
cj)i(t) for i = 1,2,3 with t and t related through equation (I19al) . There exists 
a unique maximal solution ( N,P,Z) T : [0,ic) —► R 3 , where w £ (T, oo]. Using 
the results in ©, ( N,P,Z) T : [0,wj) —> R 3 defined by P(t), Z(t)) T = 

(N(t ), P(t ), Z(t)) T is the unique maximal solution to the initial value problem 

©)-(© with w = lim t -_^_ / 0 ‘i?*(i?(P(r))) _1 dr. 

If w = oo, then 

t r r* , 

w = iim / --- dr = oo 

t—>oo J o R(P(r)) 

since R*/R{P{r)) > R*/R(N T ) > 0. 

If w < oo, then we have that lim^^- (N(t), P(t), Z(t)) T £ dCl [IT]. Since 
(N(t),P(t),Z(t)) T is bounded by Proposition 15.11 in order for the solution to 
approach dfl we must have that 

_lim R(P(t)) = 0. 

t—yw 

This implies that lim^,-- P(t) = 0. From equation (122bl) and the properties of 
h and R in equations Q and ©, there exists a positive A such that 

> -A 
at 

for t £ (w — T,w). It then follows that R(P(t)) < A(w — t) for t £ (w — T, ui), 
which implies that 


w = lim 

t—>W~ , 


f R* 

/ --- dr > iim 

'o R(P(r)) i-ew~ 


R* 


li-T A(u) - r ) 


dr = oo. 


The last part, that N(t), P(t), Z(t) £ (0 ,Nt) for all t £ [0,oo), follows di¬ 
rectly from Proposition IS.ll and the fact that ( N(t ), P(t), Z(t)) T = ( N(t ), P(t), Z(t)) T . 

□ 


6 Equilibrium Solutions 

To begin, note that (N*, P*, Z*) T is an equilibrium solution of the DDE model 
(1^1) if and only if ( N*,P*,Z*) T with R(P*) ^ 0 is an equilibrium solution 
of the TDE model (fl5l) . For a fixed value of Nt, the system (12^1) has an 


10 









equilibrium solution (TV(t), P(t), Z(t)) T = ( N*,P*,Z*) T for t > 0 if TV*, P*, Z* 
are constants that satisfy 

-^y[-/zP7(TV*) + AP* + SZ * + (1 - 7 )gZ*h(P*) 

(27a) +6 0 (N t -TV* - P* - Z *)] = 0, 

(27b) - AP* - gZ*h(P*)\ = 0, 

(27c) ^[ 75 e-W*(^Z*MP*) _ SZ *] = 0. 

If <5o = 0, then these three equations are redundant in the sense that if two 
are satisfied, then the third is also true. To avoid this problem, we use the 
conservation law m in place of equation (127al) . For P* ^ 0, equilibrium 
solutions should satisfy 

1 _ p-Som/R(P*) 

(28a) N* + P* + Z*+'ygZ*h{P*) --- N T = 0, 

o 0 

(28b) gP*f(N*) - AP* - gZ*h(P*) = 0, 

(28c) 75 e“ 5om/fl(p * ) Z*h(P*) - SZ* = 0, 

with —-—- = m/R(P*) when Sq =0. If Sq ^ 0 and P* ^ 0, equa¬ 

tions (127a| - (127cP are satisfied if and only if equations (I28al) - (l28cl) are satished. 
However, when <5o = 0, equations (I28a[) ~ (128c[) give the equilibrium solution corre¬ 
sponding to the value of total biomass, Nt, in which we are interested, whereas 
equations (127a[) - (|27c|) do not. 

There are two types of equilibrium solutions that can exist: E\ = (Aq*, P*, 0) T , 
and P 2 = (ACj, P 2 *, Z^) 1, ■ We will say that an equilibrium point exists if all its 
components are non-negative, as these represent physical quantities. We will 
only consider Z^ > 0 to distinguish between the two types. 

Note that any equilibrium solution (TV*, P*, Z*) of the DDE or TDE model 
corresponds to an equilibrium solution (TV*, P*, Z *, p*{s)) of the PDE 1 or 
PDE2 model, with 


P*(s) 


igZ*h(P*) 

R{P*) 6 


Wp* 7 . 


The PDE1 and PDE2 models admit the solution (TV, P, Z 1 p(. : s)) T = (TV^, 0, 0,0) T 
We cannot define a corresponding solution (TV, P, Z) T = (Nt, 0,0) T to the TDE 
model or the DDE model (12^1) since neither model is defined for P = 0. 
However, it is possible for these models to asymptotically approach the point 
(TVr,0,0) T . For this reason, we will define the point eg = (TVr,0, 0) T , where 
the lower case e is for emphasis that this is not an equilibrium solution, but a 
limit point. 
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6.1 The Phytoplankton-only Equilibrium E\ 

Considering Ei, we see that Z* = 0 implies that (I28cl) is satisfied. Then, for 
positive P*, (I28bl) is satisfied if and only if N* = f~ 1 (X/g). Consequently, (I28al) 
is satisfied if and only if P* = N T — N Therefore, E\ exists if and only if 

Nt > Nti, 

where 

(29) Nti = f~ l (£) . 

Note that as Nt increases, P* increases linearly while IV* remains fixed. That is, 
while Ei is a stable equilibrium, increasing the total biomass in the ecosystem 
increases the phytoplankton population while the dissolved nutrient remains 
fixed. 


6.2 The Phytoplankton-Zooplankton Equilibrium E 2 


For the equilibrium E 2 , we have that Z 2 > 0, so it is required that 

ige - 5orn/R ( p ^h(P£) -6 = 0, 


in order for (128cl) to be satisfied. This is true if and only if 


(30) 


®T) 

6 


( 7 gh(P 2 *) ^ 


Given to, P 2 is defined implicitly through this equation. Under the assumption 
that R and h are both increasing and saturating functions, we can see that 
P 2 increases with to, P 2 * = h~ 1 (6/'yg) when m = 0, and P 2 * -> oo as m Z 1 
Poo \n(^fg/8)/5 q. There is no solution for to larger than this value. However, if 
m € [0, Poo ln(7<7/5)/<So), then there is a unique positive value P 2 * that satisfies 
®. 

Assuming that to £ [0, Poo In (”/g/5)/5o), we find that (128bl) is satisfied if and 
only if 


^2 = (nWZ) - A) 


gh(P*)' 


For Z 2 > 0 , we require iV 2 > Nti- Then (I28al) is satisfied if _/V 2 satisfies 


N; + P 2 * + (gf(N*) - A) 


Pf 


gh(PZ) 


1 + 7 gh{P%) 


l _ e ~Som/R(P2) 


— Nt = 0, 


with P 2 fixed and given by ®. We can see that N£ increases with increasing 
N t and that _/V 2 = Nti when Nt = Nti + P 2 ■ Therefore, E 2 exists and is 
unique if and only if to £ [0, Poo In ('yg/S)/8 q) and Nt > Nt 2 where 


Nt2 



+ P2, 
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Parameter 

Value 


5.9 day” 1 

A 

0.017 day” 1 

9 

7 day -1 

7 

0.7 

<5 

0.17 day -1 

fW 

N 

N+k 

h{P ) 

P 

P+K 

k 

1.0/iM 

K 

1.0/xM 


Table 1: Parameter values used for all computations 


with P 2 defined through equation (1501) . Since Nt 2 > Nti, the existence of E 2 
implies the existence of E\. As Nt is increased, P 2 is fixed while iV| and Z 2 
increase. 

Figure |T] shows how the dominant equilibrium solutions change as the total 
biomass increases. We use parameter values given in Table 16.21 which were 
taken from [18]. We consider the case where R(P) = P/(P + 0.159) and various 
values of m. The top three plots are for So = 0 and the bottom three are for 
<5o = S. We see that N* = Nt and P* = Z* = 0 for Nt < Nti (eo is plotted). 
Then for Nti < Nt < Nt 2 (Pi is plotted), P* increases linearly with Nt while 
N* is fixed at Nti and Z* is fixed at 0. Then for Nt > Nt 2 (E 2 is plotted), P* 
is fixed and N* and Z* increase with Nt , though Z* saturates and N* increases 
indefinitely. Note that Nt 2 depends on m when <5o ^ 0, but not when <5o = 0. 
Also, P* is independent of m when <5o = 0, but changes for Nt > Nt 2 when 
So ^ 0. The value to which Z* saturates does not change with m when So = 0, 
but increases with m when Sq ^ 0- 


7 Linearization and Characteristic Equation 


In this section, we will show that a general threshold delay equation and its 
corresponding delay differential equation, as in |lf)j , have identical linearizations 
about an equilibrium solution. 

Consider the threshold delay equation in the general form 

( 31 ) —jp- = f{x{t),x(t-T(xt)),T(x t )), [ K(x t (u))du = m, 

dt J-t(x t ) 

and the corresponding functional differential equation 


( 3 2 ) 

““TT = x(t - m/K *), t(%)), 

at K(x(t)) 


Hxf) 



K * 

K(xi(u)) 


du , 


where / and K are differentiable, and K* is a typical value of K. Results on 
the linearization of state-dependent delay differential equations with bounded 
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Figure 1: Equilibrium solutions as a function of total biomass for m = 
0,5,10,15,19.7 and R(P) = P/{P + 0.159). The left three are for S 0 = 0 
and the right three are for 5o = <5- If £2 exists, it is the solution plotted. If E 2 
does not exist, but E\ does, then E\ is plotted. If neither exist, the limit point 
( Nt , 0,0) T is plotted 
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delay can be found in [J. These results apply to systems with bounded delay, 
which do not apply necessarily apply directly to the system here. Nevertheless, 
we will show that applying the linearization procedure for a state-dependent 
delay system with bounded delay to system (1311) yields the same linearization 
as linearizing system m in the standard way for fixed delay systems. 

Let x* be an equilibrium solution of equations (13T1) and (13^1) . That is, a 
constant function such that f(x*,x*,r*) = 0 where r* = m/K(x*) = t(x*) = 
t(x*). Let y(t) = x(t) — x*. By setting the state-dependent delay to its equi¬ 
librium value, and linearizing equation (1311) as in [3], we obtain 

= D d{x*,x*,T*)y{t) + D 2 f(x*,x* 1 T*)y(t-T*) 

(33) + D 3 f(x*,x*,T*)Dr(x*)y t 


where Dif is the derivative of / with respect to its ith argument and 

(34) Dr(x*)yt = /° T , yt{u)du - 

Let y(t) = x(t) — x*. Then linearizing equation (1321) in the standard way 


yields 


(35) 


dy(t) 

dt 


= T -( H, A D if( x *, x *, T *)y(i) + D 2 f{x*,x*,T*)y(i-m/K*) 
K [x ) 

+ D 3 f(x *, x* , t*)Dt(x* )y t -] 


where 

(36) Df(x*)y { = ~ K Kfrxt ^ [ Vt( u ) du - 

^ [X ) J—m/K* 


Without any loss of generality, we can set K* = K(x*). Indeed, this would 
be the natural choice of K* if we were interested in solutions near a given 
equilibrium. Then it can be see that equations (1331) and (1351) are identical. 

We can check that (l34l) is correct by verifying that 


(37) 


lim ■ 

fc.-> o 


t( x* + h) — t(x*) + 


DK(x*) fO 
K(x») 


f _ T » h[u) du 


= 0 , 


where h £ C[— r, 0] for some r > for all <f> in a neighbourhood of x*. Here, 
||.|| is the usual sup norm. We proceed as follows: 


t(x* + h) — t(x*) + 


DK(x* 


h(u) du 


K(x*) 


K(x*) 

/ 0 rO pO 

K(x*)du— / K(x*) du + DK(x*) / h(u) 

-r(x*-\-h) J—r* J—r* 


du 
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Under the assumption that K £ C 1 , it is true that K{x*) = K{x* 
DK(x*)x + G(x), with G satisfying 

lim yprp 1 = 0, 

x ~>° |N N 

where ||.||e is the Euclidean norm on R". Then 


+ x) - 


' -r(x*-\-h) 
r 0 


K(x*) du 


[ K(x* + h(u)) du — DK(x*) [ 

J-r(x*-\-h) J — 


' — t{x*-\-H) 

+ f G(h(u)) du. 

J—t(x* -\-h) 


r{x*-\-h) 


h(u ) du 


Then by the definition of r: 

r0 

r(x*-\-h ) 


/ 0 pO 

K(x* + h(u)) du = / K(x*)du = m. 

-t(x* -\-h) J —r* 


So it can then be 

seen that 

lim 

h^>- 0 

t{x* + h) - 


= lim 
h^-0 

DK(x*)fZ 



DK(x*) r 0 


DK(x*) 5 T T f +h) h(u) du + f" T(x , +h) G(h(u)) du 


K (ad 


= 0 , 


from the properties of G. 

To verify equation (Idfil) . we check that 


(38) 


lim 
h —>-0 


t(x* + h) — t{x*) 


K* DK(x*) r 0 
K(x *) 2 J—m/K* 


h{u ) du 


= 0. 


This follows from the definition of r in equation (1321) and the fact that for x 
near x* 


(39) 


1 1 DK{X *K I C(x) 

K(x* + x) K{x*) K{x*) 2 { h 


for some function G that satisfies 


]G(®)I 

INI 
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In particular, the linearized system for our model l (fl5l) or (021) 1 about an 
equilibrium solution (TV**, P*, Z*) T with R* taken to be R(P*) is 


(40) 

where 

A 1 = 


dy(t) 

dt 


= A\y(t) + A 2 y(t — T) + A 3 j y(t + u)du 


—yP*a — — yc + A + (1 — 7 )gZ*b — 60 6 — So + (1 — 7 )gd 


yP*a 

0 


/re — A — gZ*b 

e-^ igZ *d«!f 


-gd 

-6 


Ao — 


0 e ~ s ° T 1 gZ*{b-^fdd) e~ 5 ° Tr ygd 


A, = I 0 


0 0 

0 0 

-S n T„,„7*jR'{P') 


0 S 0 e- d oT igZ *d^f 0 

and a = /'(IV*), b = h'(P*),c = /(IV*), d = h(P*), and T = m/R(P*) = 
t(tti, P*) = r(m, P*). 

Substituting y{t) = ve st into (OU1) . with u £ l 3 , we can obtain the charac¬ 
teristic equation: 


(41) 


det ( si — Ai — A 2 e s + A, 


(l-e- sT ) 


= 0 . 


If all values of s that satisfy equation m have a negative real part, then 
(. N* 1 P*,Z*) t is an asymptotically stable equilibrium solution of system (1221) 
m- Due to the equivalence of systems m and (l22l) . they are asymptotically 
stable under identical conditions. Thus, we can use the characteristic equation 
m to obtain asymptotic stability in system (ED- Therefore, if all values of s 
that satisfy equation (PHI) have a negative real part, then (IV*, P* , Z*) T is an 
asymptotically stable equilibrium solution of system (ED- 


8 Stability and Extinction 

8.1 Conditions for the Extinction of Plankton 

It is possible for the solutions to the TDE model (fl5l) to approach the point 
eo = (Nt, 0, 0) T . To show this, we will use the following Lemma. 

Lemma 8.1 If (N(t ), P(t ), Z(t)) T is a solution to system (1151) such that 
lirn^oo P(t) = 0, then lim t _> 0 0 (N(t), P(t), Z(t)) T = (N T , 0, 0) T . 


17 







Proof: Assume that lim^oo P(t) = 0. The term e _l5 ° T ( m ’ Pt ) is clearly 

bounded. Also, by Proposition 15.21 Z(t — T(m,P t )) < Ar, so this term is 
also bounded. From the properties of h and R in equations (U) and <[5|), there is 
a function, w, such that lim^oo w(t) = 0 and dZ J t t ' > < w(t) — 5Z(t) for t > to- 
integrating this inequality, we obtain 

Z{t) < e - 4 (*-*o>Z(to) + e~ st f e Su w{u) du. 

■'to 

In the case where e Su w(u) du < oo it is obvious the final term approaches 
zero as t — > oo. Otherwise, we can apply l’Hopital’s rule and get the same result. 
Since Z(t) >0 by Proposition 15.21 we have that lim^oo Z(t) = 0. 

The result then follows from the conservation law m if it is true that 

lim r ^ = „ 

t^ooj o R(P(t-T(s,P t ))) 

If lim^oo r(s, Pt) < oo, then Z(t — r{s,Pt)) —>■ 0 and the result immediately 
follows. If lim^oo r(s, Pt) = oo then e~ s ° T ^ s,Pt ^ —> 0 and the result follows. □ 
This also implies that if P(t) in the solution to the DDE model (l22l) reaches 
zero in finite time, then Z(t ) reaches zero at the same finite moment in time. 
This is because the transformation maps t = oo to t = w< oo. 

The following Proposition shows conditions for the extinction of plankton. 

Proposition 8.2 If 0 < Nt < An, where An is given in (l29l) . and fa, <fz) T € 
Dn t , then the solution to the initial value problem (EMU asymptotically ap¬ 
proaches (Ar,0,0) T . 

Proof: From (Il5bl) . we have that < pf {N^Pit) — \P{t). Then Nt < An 
implies that /z/(Ar) — A < 0, which implies that P{t) < P(fo)e ( - M ^ JVT ) _A )( t_ *°). 
Since Pit) > 0 from Proposition 15.21 lim t ^. OCJ P(t) = 0. The result then follows 
from Lemma ED □ 

8.2 Global Stability of E { 

Here we will show conditions for the global stability of . Consider the follow¬ 
ing Lemma. 

Lemma 8.3 Let (</>i, 4>3) T G Dn t and denote ( N,P,Z) T : [0,oo) —> K 3 as 
the solution to the initial value problem CHD-CED. If Nt > An, then for any 
constant /3 > Nt — An, there exists at\ >to such that P(t) < j3 for all t>t\. 

Proof: Fix f3 < Nt — An- From the conservation law m and Proposition 
15.21 A(f) < Nt — P(t). Then for P(t) > j3 we have 

^<W(A(t))-A]P(t), 

< \pf (N t - P{t)) - A]P(t), 

< [pf[N T -P)- A }P(t). 
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Since h/(Nt — /3) — A < nf(Nxi) — A = 0, we have 

^ < \pf(N T - /?) - A]/? 

for P(t) > (3. Therefore there exists t\ > to such that P(t) < (3 for all t > ti. □ 

Proposition 8.4 Let </> 2 , </> 3 ) T £ Dn t and denote ( N , P, Z) T : [0, oo) —> M 3 
as the solution to the initial value problem (1141) - (1151) . If Nt i < Nt < Nt 2 , then 
lim i ^ 00 {N(t) 1 P(t) 1 Z(t)) T = (Nf,Pf, 0) T . 


Proof: Note that /if(Nf) — A = 0 and P* = Nt — N£ < P 2 * for Nt < Nt 2 - 
Consider the function 

r N ( t '> r p W x + m - N t 

V\{t) = / [/j,f(x) — \]dx + 5o / - dx. 

Jn * Jn t -n * * 

We have that V\ (t) > 0 for P(t) > 0. Taking the derivative along solutions to 
the TDE model (IT31) . we obtain 

jVi(t) =[pf(N(t)) - A ][-gP(t)f{N{t)) + A P(t) + 8Z(t) + (1 - 7 )gZ(t)h{P(t)) 
+ 5 0 (N T -N(t)-P{t)-Z(t))] 

+ So m+ p^~ NT [tiP(t)f(N(t)) - A P(t) - gZ(t)h(Pm, 

< - P{t)[pLf(N(t)) - A ] 2 - S 0 [N(t) - Nf}fif(N(t)) - A] + MZ(t), 


for some positive constant M. This M exists because N(t) and P(t) are bounded 
and from the properties of h in (U)- Then consider the function 


V 2 {t) 


Z(t) + 


f 

J t—T(m,Pt) 


7 ge- 50 ^™'^Z(u)h(P(u)) du. 


We have that V 2 (t) > 0. Choose f3 £ ( Nt — Nti,P 2 ). Then by Lemma IQ1 
there exists t\ > to such that P(t) < [3 for t > t\. This implies that 

7 ge- SoT ( m ’ Pt) h(P{t)) < ^ge- s ° T{m ’^h{f3) < 7 ge- s ° T ^ m ’ P ^ h(P^) = S 


for t>t\. Define the positive constant a = S — 7 ge SoT ( m ’^h(/3). Differentiat¬ 
ing, we have that 

- «<*> 

+ 7 ge- SoT ^ Z(t)h(P{t)) 

- 7 ge-^^Zit - r(m, P t )h(P(t - r(m, P t )) (l - , 

< — aZ(t), 
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where we have used that 


l dr _ R(P(t)) 

dt R{P{t — T(m, Pt)) ’ 

h{P(t)) < h(/3), and e~ SoT ( m ^ < e -^T(m,/3) for t > tl s etting y{t) = Vi (t) + 
M±Ml y 2 where Mi > 0, we have that V(t) >0 and that 

jV(t) < -S 0 [N(t) - NfiW(N(t)) - A] - M l Z{t), 

for t, > t, i. Integrating, we obtain that 

(42) V(i) + [ [5 0 (N(u) - - A) + M^u)] du < V(h). 


Since V (t) and the above integrand are nonnegative, and since the inequality 
(l42l) must be true for all time, it follows that 


'to 


[5q(N(u) — Ni)(/j,f(N(u)) — A) + MiZ(u)] du < oo. 


By Barbalat’s Lemma 0, it must be true that N(t) —> N* and Z(t) —> 0. By 
the conservation law (TT51) it then follows that P(t) —> P£. □ 

The attractivity of E\ has thus been shown. However, Proposition 18.41 does 
not address its stability. Let ( N , P, Z) T : [0, oo) —► R 3 denote the solution to the 
initial value problem (fl4l) - (fl5l) . We will consider an equilibrium point E to be 
asymptotically stable if for any e > 0, there exists a neighbourhood N e around 
E such that if (^ 1 ,^ 2 , ^ 3 ) £ P>n t C\ N e then ||(7V(t), P(t), Z(t)) T — E || < e for 
t > t 0 and limt->oo{N(t), P(t), Z(t)) T = E. 

We will use the linearization of the TDE model (fl5l) to prove the following 
Proposition. 


Proposition 8.5 If Nti < Nt < Nt 2 , then E\ is asymptotically stable. 


Proof: The condition Nti < Nt ensures that E\ = (1V*,P 1 *,0) T exists. The 
characteristic equation of the corresponding linear system is 


( s + pPfa + S 0 S 0 
—pPfa s 

0 0 

= (s + [iPfa)(s + Sq){s + S 


-S + S 


0 - (1-7 )gd 

gd 


s + S — 'ygde (Ao +s)t 
7 gde- {5o+s)T ) = 0 . 


The first two factors give two negative real eigenvalues. Since N T < N T2 , it is 
true that r )gde~ &aT < S. Set s = a + iuj for real a and ui and assume that a > 0. 
We have that s+<5— r )gde~ <K&0+s ' ,T = 0 if and only if a+iui+8 = 7 gde~^ 0+a+lul ' >T . 
Squaring the modulus of both sides we obtain 

(a + S) 2 +uj 2 = { 1 gd) 2 e- 2{5o+a)T < S 2 , 
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which implies a < 0. This contradiction shows that there are no solutions with 
nonnegative real part. Since all the eigenvalues have negative real parts, E\ is 
asymptotically stable. □ The global 

attractivity together with the local stability of E\ for Nti < Nt < Nt 2 shows 
that Ei is globally asymptotically stable on Dn t . 

8.3 Stability of E 2 : No Delay 

If m = 0 then the zooplankton are considered mature immediately at birth. In 
this case there is no delay and system m is just a system of ordinary differential 
equations. In El, we showed that if h has a negative second derivative (a Type 
II response) and / has a negative second derivative, which is the case we will 
consider in the following sections, there is a unique Nt 3 > Nt 2 such that E 2 
is asymptotically stable if Nt 2 < Nt < Nt 3 and unstable if Nt > Nt 3 - This 
value of Nt 3 is independent of Sq and the functional form of R{P) , since m = 0 
implies that Sq and R(P) do not play a role in system (1151) under the assumption 
that c/> 2 , fa) € Dn t . 

8.4 Stability of E 2 : State-Independent Delay 

Consider the case where R(P) is constant. That is, the growth rate of zooplank¬ 
ton is independent of the phytoplankton population. Without loss of generality, 
we will assume m has been scaled so that R(P) = 1. In this sense, m has dimen¬ 
sions of time, and represents the age at which zooplankton reaches maturity. In 
this case, the delay is no longer state-dependent and is fixed at r = m. 

Due to the complicated nature of the resulting characteristic equation, we 
use numerical methods to study the stability of E 2 (both here in the state- 
independent delay case and in the following section, which deals with the state- 
dependent delay problem). We use the parameter values in Table 171721 

By setting s = ico in the characteristic equation m and allowing the m 
and Nt parameters to vary, we can find solutions where the eigenvalues have 
zero real part. These represent critical points where stability might switch. 
Since the characteristic equation is complex valued, it gives us two equations 
that need to be satisfied. The equilibrium equations change with changing 
m and Nt, so (l28l) gives us three more equations to solve. We then have 
six unknowns: N* , P* , Z* ,m, Nt,w. We can then find one-dimensional curves 
where the five equations are satisfied, using pseudo-arclength continuation [9]. 
Note that if (N* , P* , Z* , to, Nt,uj) solves the five equations then w represents 
the frequency for solutions to system m near the corresponding equilibrium 
solution (N* , P* , Z*) when the size at which the immature zooplankton reach 
maturity is m and the total biomass is Nt- 

Figure [2] shows regions in the m — Nt plane where the equilibrium solutions 
exhibit different behaviour. This was done for two cases: Sq = 0 (left in Figure 
0 and Sq = S (right in Figure 0. The former case represents the situation 
where the immature zooplankton have a zero death rate. That is, we ignore the 
possibility of them dying before reaching maturity. In region 1 we have that E\ 
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Figure 2: The solid curves are where the linearized system has an eigenvalue 
with zero real part for R(P) = 1 and do = 0 (left) and do = d (right). Region 
1 is where Ei and E 2 do not exist. Region 2 is where Ei exists, but E 2 does 
not. Region 3 is where E\ and E 2 exist. The subset of region 3 under the solid 
curves is where E 2 is asymptotically stable 


and E 2 do not exist. This region is the same, regardless of the value of do- In 
region 2, E\ exists, but E 2 does not. This region does depend on the value of 
do, as a larger do requires more biomass for the E 2 equilibrium to exist. Region 
3 shows where Ei and E 2 exist. We note that Roo In( 7 p/d)/do ~ 19.77, which 
is an upper limit for m when do > 0. The solid curves shows the values of m 
and Nt where the linearized system has an eigenvalue with zero real part. The 
subset of region 3 below the solid curves represents values of m and Nt where 
the E 2 equilibrium is asymptotically stable. This follows from the continuity 
of the eigenvalues with respect to parameter values and the fact that there are 
no eigenvalues with nonnegative real parts when m = 0 and Nti < Nt < Nts 
(between where the upper dotted curve intersects the vertical axis and where 
the solid curve intersects the vertical axis). 

These regions are all fairly independent of m when do = 0. That is, if the 
juvenile zooplankton have a zero death rate, then the required level of maturity 
has little effect on the stability of the equilibrium solutions. However, when 
there is positive death rate for the immature zooplankton, the required level 
of maturity plays a more important role. We see that critical values of total 
biomass increase with increasing m. 

Figure [3] shows the frequencies that correspond to the curves in Figure [2j 
These are the imaginary parts of the eigenvalues with zero real parts along 
the curves. For values of m and Nt near the solid curves in Figure [2j we 
would expect a slow growth or decay rate in the solution to system m when 
it is near the equilibrium solution, and for it to have a frequency close to the 
corresponding value of ui in Figure [3j 
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l = 0,S o = 0 


l = 0, Sq = 5 



Figure 3: The corresponding frequencies to the curves in Figure [2] The value of 
w is the imaginary part of the eigenvalue with zero real part along these curves 


8.5 Stability of E 2 : State-Dependent Delay 

Here we consider the case where R(P) = -prj for various values of l. Using the 
same numerical technique as in the previous section, we compute curves in the 
m — Nt plane where there is an eigenvalue with zero real part. This tells us 
parameter values where stability can change. 

Figure |4] shows regions in the m — Nt plane that have different behaviour. 
This was done for <5 0 = 0 (left figures), Sq = S (right figures) and for l = 0.01 
(top figures), l = 0.159 (figures second from the top), l = 1.00 (bottom figures). 
As before, region 1 is where neither Ei nor E 2 exist. Region 2 is where Ei exists, 
but E -2 does not. Region 3 is where both E\ and E 2 exist. The solid curves are 
where there is an eigenvalue with zero real part for the system linearized about 
E 2 . Hence, the subset of region 3 beneath the solid curves represents the values 
of m and Nt where the equilibrium solution E 2 is asymptotically stable. 

We can see that the case where l = 0.01 is very similar to the case R{P) = 1, 
shown in Figure [2] This is due to the fact that l in this case is small enough 
relative to P so that R(P) is approximately constant. 

When l is increased to 0.159, its value is no longer small relative to typical 
values of P (see Figure [TJ. In this case, the delay has a stronger dependence on 
the quantity of the phytoplankton. 

We see many more curves in the case of <5o = 0. The minimum values of Nt 
for each of these curves increases slightly as m is increased. For l = 1.00 there 
are about an order of magnitude more curves. The minimum values of Nt for 
these curves increases more as m is increased, which creates an overall larger 
region of stability than when l = 0.159. 

For l = 0.159 and So = 6 there is a single curve that loops three times. In 
fact, we observed that l = 0.01 is a single curve with three loops as well, but 
we would have to plot it for Nt values much larger than 100 in order to see 
this. As we vary l from 0.01 to 0.159, we see this loop become tighter. As l 
increases to 1.00, the loops either vanish or become smaller than what can be 


23 




detected at this scale. Either way, we get a much simpler region of stability, 
even though a larger value of l means that the delay has a stronger dependence 
on the quantity of phytoplankton. 

Figure [5] shows the corresponding frequencies to the curves in Figure 01 In 
this figure, w is the imaginary part of the eigenvalues with zero real parts. In 
the case where So = 0, it is true that the characteristic equation is periodic in 
to in the sense that if the characteristic equation is zero when m = mo, then 
it is also zero when m = mo + 2mrR(P*)/oj. This is because m only appears 
in the characteristic equation in the form e imu / R ( p ). This periodicity can be 
easily seen in Figure [5] when <5o = 0. Since changing l by an order of magnitude 
changes R(P*) by an order of magnitude, we get an order of magnitude more 
curves for the same range of m. 


9 Numerical Simulations 

In order to verify some of the results from the previous section, we performed a 
series of numerical simulations. Rather than simulating system (usd directly, we 
simulated (l22l) and then transformed the resulting time series using equations 
(U9bl) - (ll9dl) and (12T1) . System (l22l) was simulated with a second-order method, 
which included numerically solving the integral (l20l) using a trapezoidal rule. 

The simulations allowed us to verify the regions of stability for E 2 that were 
computed numerically. For instance, we varied values of m and AY near the 
solid curves in Figure 0J The initial conditions for the simulation were chosen 
close to the equilibrium solution E 2 . For values of to and A'Y that were in a 
stable region, we verified that the time series of the simulation decayed in time. 
For values of to and AY that were in an unstable region, we verified that the 
time series of the simulation grew in time. Many tests were done for various 
values of l and for <5o = 0 and <5o = 5, and no inconsistencies were found. That 
is, the simulations always agreed with the numerical bifurcation analysis. 

Figure [6] shows one such verification. For l = 0.159, So = S , and to = 6, 
Figure 0] predicts that E 2 should be stable for AY = 10 049 and unstable for 
AY = 10° 51 . On the left in Figure 01 the simulation for AY = 10° 49 indeed 
suggests that the equilibrium solution is stable. On the right in Figure 01 the 
simulation for AY = 10°’ 51 suggests that the equilibrium is unstable. While 
such simulations are by no means a proof of stability, it is still reassuring that 
they agree with the predictions from the numerical bifurcation analysis. 

Figure 0] shows a subinterval of the time series of a simulation for l = 0.159, 
So = 5, to = 8, and AY = 10 0,73 . For these values Figures 0] and 0] together 
show that there are about five pairs of eigenvalues with real parts close to zero. 
Thus, we might expect the solution to exhibit up to five frequencies. In Figure 
0] we see that the solution is very irregular, which is not surprising given the 
predicted spectrum shown in Figure 0] 
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I = 0.01, do = o 


l = 0.01,5 0 = J 



Figure 4: The solid curves are where the linearized system has an eigenvalue 
with zero real part for R(P) = ppj for various values of l and S 0 = 0 (left) and 
do = 6 (right). Region 1 is where E\ and E 2 do not exist. Region 2 is where 
Ei exists, but E 2 does not. Region 3 is where E\ and E 2 exist. The subset of 
region 3 under the solid curves is where E 2 is asymptotically stable 
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i = O.Ol,0o = O 


l = 0.01,5o = 0 



Figure 5: The corresponding frequencies to the curves in Figure [2] The value of 
u) is the imaginary part of the eigenvalue with zero real part along these curves 
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Figure 6: The time series of a simulation for R = P+9 159 , 5o = 6, m = 6, 
and Nt = 10°' 49 (left) and Nt = 10° 51 (right). The initial conditions were 
chosen near the equilibrium solution. In the case of Nt = 10° 49 , the simulation 
suggests that the equilibrium solution is stable, but in the case of Nt = 10 a51 , it 
suggests that the equilibrium solution is unstable. This agrees with Figure[4](on 
the right, second from the top), which indicated that there is a zero eigenvalue 
at to = 6 and Nt = 10 0,5 °. Also plotted is the time series of the r(m, Pr), 
which is the state-dependent delay 



Figure 7: A subinterval of the time series of a simulation for R = P+9 159 , 
<5o = S, to = 8, and Nt = 10 0 73 . The state-dependent delay, t(to, Pt) is also 
shown. This simulation shows the possibility of very irregular and possibly 
chaotic solutions to system da 
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10 Discussion 


We have looked at a model of a closed planktonic ecosystem that depends on 
the maturity structure of the immature zooplankton. Using techniques in m 
and |20j , we were able to transform the model into a delay differential equation 
with a state-dependent threshold-type delay or a delay differential equation 
with a state-independent delay. Such transformations allowed us to use results 
readily available from the theory of delay differential equations to study the 
qualitative features of the model, such as existence and uniqueness of solutions, 
boundedness, persistence, and stability. For instance, we showed that solutions 
exist for all time, remain positive, and are bounded, which are desirable features 
for an ecological model to have. 

Being able to represent a system in a variety of ways is very beneficial from 
a mathematical perspective, as we are able to choose the framework that is 
most convenient for the situation at hand. For instance, the state-independent 
delay differential equation was useful for applying well-established theory for 
functional differential equations with fixed delay, as well as performing numerical 
simulations. The state-dependent threshold differential equation was convenient 
for computing the linearization and for arguing qualitative results. The partial 
differential equation was useful since it represents the model in its most intuitive 
form, and therefore offers the best framework for interpreting results. 

Furthermore, in a practical sense, the different models offer flexibility in 
what initial data is required for numerical simulations. For instance, in some 
situations it may be convenient to measure the spectrum of the juvenile zoo¬ 
plankton at a given time, while in other situations it may be more convenient 
to measure the histories of the phytoplankton and mature zooplankton over a 
sufficiently long time interval. Since it is possible to transform between the PDE 
and TDE models, the type of data obtained for the initial conditions does not 
necessarily dictate which equations we need to use to simulate the ecosystem. 

A key parameter in the study of closed ecosystems is the amount of biomass, 
which is fixed by the initial conditions. General results include the existence 
of two critical values of the total biomass. The first is the minimum amount 
needed to sustain the phytoplankton population, which we have called Nti- If 
the biomass is less than Nti, then the phytoplankton and zooplankton both 
become extinct, as there is not enough biomass to sustain their populations. 
Conversely, if the biomass is greater than Nti, then the phytoplankton do not 
become extinct. The second critical value of the total biomass, which we have 
called Nt 2 , is the minimum amount needed to sustain the zooplankton popu¬ 
lation. If the total biomass is less than Nt 2 then the zooplankton population 
becomes extinct, but does not if the biomass is greater than Nt 2 - We have 
shown that if the total biomass is greater than Nti, but less than Nt 2 , then 
the system globally approaches a unique, phytoplankton-only equilibrium solu¬ 
tion that depends on the total biomass. Future work might include further study 
of the system when the biomass is greater than Nt 2 , including a formal study 
of weak/strong (uniform) persistence [3] and possibly global behaviour of solu¬ 
tions. Here, for Nt > Ntz, we mainly focused on local stability of the unique 
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equilibrium solution with numerical techniques for a chosen set of parameter 
values. 

In the case when the juvenile zooplankton have a zero mortality rate, (So = 
0), the stability results have a strong dependence on the total biomass (Nr) 
while the required level of maturity (m) is less significant. In essence, it seems 
that the time to maturity does not matter if juvenile zooplankton are not being 
lost due to mortality. However, when we allow them to have a positive mortality 
rate (<5o > 0), the required level of maturity becomes more important with 
regards to stability. In this case, Figure 0] suggests that m must increase with 
Nt in order to maintain stability of the phytoplankton-zooplankton equilibrium. 

As well as stable behaviour, the model ecosystem can simple exhibit periodic 
solutions as well as more complicated dynamics. We saw complicated orbits in 
many cases where parameters were such that the characteristic equation had 
multiple pairs of eigenvalues that had real parts close to zero. 

Further work may include adding structure to the immature phytoplankton 
population in a similar way to what was done for the immature zooplankton. In 
|13] , we have included a delay in nutrient recycling while ignoring the structure 
of the immature zooplankton, but a more complete model would include both 
of these effects. It may also be worthwhile to see how non-linear closure terms 
affect the overall behaviour of the system when the zooplankton size-structure 
is present. An investigation into the role of this closure term in models without 
such size-structure is given in [5] and is generally considered to be very significant 
in determining the dynamics that can occur. Spatial structure can also be 
added, but similar transformations between partial differential equations and 
delay differential equations may not be possible in this case. 
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